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In the present work, we investigate how single- and multi-vortex-ring states can emerge from a 
planar dark soliton in three-dimensional (3D) Bose-Einstein condensates (confined in isotropic or 
anisotropic traps) through bifurcations. We characterize such bifurcations quantitatively using a 
Galerkin-type approach, and find good qualitative and quantitative agreement with our Bogoliubov- 
de Gennes (BdG) analysis. We also systematically characterize the BdG spectrum of the dark 
solitons, using perturbation theory, and obtain a quantitative match with our 3D BdG numerical 
calculations. We then turn our attention to the emergence of single- and multi-vortex-ring states. 

We systematically capture these as stationary states of the system and quantify their BdG spectra 
numerically. We find that although the vortex ring may be unstable when bifurcating, its instabilities 
weaken and may even eventually disappear, for sufficiently large chemical potentials and suitable 
trap settings. For instance, we demonstrate the stability of the vortex ring for an isotropic trap in 
the large chemical potential regime. 


I. INTRODUCTION 


Over the last two decades, there has been an intense in¬ 
terest in nonlinear matter waves in the context of atomic 
Bose-Einstein condensates (BEGs) [H-Q. In the three- 
dimensional (3D) setting, arguably, one of the most pro¬ 
totypical excitations that arise are vortex rings (VRs). 
Such toroidal-shaped vortices, or rings of vorticity, have 
been predicted theoretically and observed experimentally 
in numerous studies; see, e.g., the relevant reviews [j- 
In addition to their relevance in superfluids, VRs 
are of intense interest in other areas, e.g., in fluid me¬ 
chanics 00. Superfluid VRs have been observed ex¬ 
perimentally in helium , far before their emergence 

in atomic BECs of dilute alkali gases. In the atomic 
BEG setting, these states have arisen in a variety of 
ways. One such example is through the bending of a 
3D vortex line (i.e., a line of vorticity) so that it closes 
on itself [HI- Vortex rings have also been experimen¬ 
tally realized through the decay of planar dark solitons 
in two-component BEGs [T^ , by density engineering^^ 
(related to earlier theoretical proposals of Refs, (l^ IlSll) . 
and in the evolution of colliding symmetric defects |I6l |. 
They have also been detected through their unusual col- 
lisional outcomes of structures that may appear as dark 
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solitons in cigar-shaped traps 0. 

Numerical studies have also theoretically explored po¬ 
tential methods to generate VRs. Some of these exam¬ 
ples involve the flow past an obstacle [H, [l^, Bloch 
oscillations in an optical trap [l^, the collapse of bub¬ 
bles [2l| , the instability of two-dimensional f^D /rarefac¬ 
tion pulses [i^ , the flow past a positive ion [2^ [2^ or an 
electron bubble [2^, the Grow instability of two vortex 
pairs [ 2 ^ and the collisions of multiple BECs VR 

cores have an intrinsic velocity |29| (a feature that can 
be understood by considering their cross-section resem¬ 
blance to a vortex dipole which is well-known to travel at 
a constant speed 0 , 0 ). However, this intrinsic velocity 
can be counterbalanced by the presence of an external 
trap, which is commonplace in atomic BECs [s^. This, 
in turn, creates the possibility for the existence of station¬ 
ary VRs, and, in facE even of multi-rings, as discussed 
earlier, e.g., in Ref. [^ . 

The aim of the present work is to revisit the study of 
VRs, but from a different perspective to the above works, 
namely we focus of their bifurcation and emergence in the 
weakly nonlinear limit. In this sense, our study is par¬ 
tially connected to the recent work of Ref. , where a 
somewhat similar approach was adopted. Nevertheless, 
their work focused on a setting where the trapping was 
only in two out of three spatial directions; another fun¬ 
damental difference with respect to Ref. is their ab¬ 
sence of stability information in the relevant setup. On 
the other hand, our study is also inspired by the pio¬ 
neering work of Ref. , which, in turn, motivated the 
experimental results of Ref. [l2|- In Ref. [s^, the spec- 
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trum of a dark soliton was examined (although not in 
its full detail, as we will discuss below) and some of the 
key instability mechanisms associated with it were eluci¬ 
dated. We are also motivated by the analysis presented 
in Refs. [13,[s^, which used a somewhat similar approach 
to that used herein to explore bifurcations from a two- 
dimensional dark soliton to vortex-dipoles and, more gen¬ 
erally, multi-vortex structures. 

Our fundamental starting point will be the realiza¬ 
tion that a VR emerges from a suitable combination of 
two distinct states with a relative phase (of 7r/2 between 
them), namely a planar dark soliton (say, along the z- 
direction) and of a ring dark soliton (RDS) along the 
(x, y) plane, another state studied extensively in earlier 
work; see, e.g., Refs, Should the planar dark 

soliton have a lower energy than the RDS, then the bi¬ 
furcation happens from the former; if the order of en¬ 
ergies is reversed, then the bifurcation occurs from the 
RDS. If, however, the energies are equal (which, as we 
will see, happens for anisotropic traps, i.e., when the 
trapping frequencies along the z- and radial r- directions 
satisfy = 20;^), the VR state already emerges at the 
linear (noninteracting) limit. The presence of two princi¬ 
pal modes in the coherent structure enables us to adapt 
a Galerkin-type, two-mode methodology, originally de¬ 
veloped for double-well potentials; see, e.g., relevant de¬ 
tails in Ref. [s^, based on the fundamental earlier work 
of Ref. [ 43 . This two-mode approach allows us to pro¬ 
vide a quantitative estimate, which turns out to not only 
qualitatively but also quantitatively predict the bifurca¬ 
tion/emergence point of the single VR. We use a simi¬ 
lar approach to generalize to more complex states such 
as the double VR (arising due to the combination of a 
double-soliton/second excited state along the z-axis, with 
a RDS along the (a;, y) plane), and we further extend this 
approach to single and double vortex lines. 

Since the planar dark soliton may spawn the VR, in the 
appropriate regime, we expand upon the work of Ref. 
by capturing all its modes using the Bogoliubov-de 
Gennes (BdG) analysis near the linear (small-interaction) 
limit. This limit enables a perturbative treatment of 
the relevant modes, both in the isotropic and in the 
anisotropic case, clearly allowing us to identify neutral 
modes, modes responsible for instabilities and also the 
bifurcations of novel states such as the VR, and the sin¬ 
gle and double solitonic vortex. We provide, both, an¬ 
alytical estimates and full-numerical calculations for the 
relevant eigenvalues, and provide connections with the 
analytical two mode theory discussed above. Once the 
relevant bifurcations arise, especially those leading to the 
single and double VRs of focal interest herein, we turn 
our attention to their spectra to quantify instability and 
the potential stabilization mechanisms. This leads us to 
conclude that these states tend to be weakly unstable in 
the small-interaction limit, but may —in principle— be 
stabilized in the large-interaction setting. 

It is relevant to mention that this problem remains 
of particular interest, not only for theoretical but also 


for recent experimental studies. A notable example is 
the very recent work of Ref. [dlj , where the tomographic 
imaging of a superfluid Fermi gas of ®Li atoms enabled 
the observation of the snaking instability of the planar 
dark soliton into a VR (and subsequently into a vortex 
line). Gonnecting state-of-the-art observations such as 
these to the understanding of stationary-state stability 
offered here should furnish a complete picture of the un¬ 
derlying complex dynamics of the system. 

Our presentation is structured as follows. Section [II] 
contains the different aspects of our analytical results; 
first, the bifurcation analysis (of single and double VRs) 
using the Galerkin method; and second, the spectral 
analysis, for the small-interaction limit, of the dark soli- 
tons which we connect to the above bifurcation structure. 
In Sec. uni we provide detailed numerical computations 
of the existence and stability of dark solitons, and sub¬ 
sequently VRs and multi-VRs, starting from the small- 
interaction limit. We provide comparisons between our 
numerical and analytical predictions. Finally, in Sec. IIVI 
we summarize our findings and present both, our con¬ 
clusions, as well as some interesting directions for future 
studies. 


II. ANALYTICAL CONSIDERATIONS 

In this work we utilize the 3D Gross-Pitaevskii equa¬ 
tion (GPE), expressed in the following dimensionless 
form [3|: 

iut =—-;^‘^u + V{r)u+\u\'^u. (1) 

Here, u is the macroscopic wavefunction of the 3D BEG 
near zero temperature, while the potential assumes the 
prototypical parabolic form, namely: 

V{x,y,z) = \u:y + \u:lz\ ( 2 ) 

The parameters iOr and Uz represent the trapping 
strengths along the (a;, y) plane and z direction, respec¬ 
tively, with the spherically symmetric (isotropic) case 
corresponding to ujr = 

A. Vortex Ring Bifurcation Near the Linear Limit 

In the linear limit, Eq. (HD reduces to the quantum 
harmonic oscillator (QHO), with energy spectrum: 

En,m,k = LOr {n + m + 1) + LOz {k + 1/2) , (3) 

where n, m and k are the non-negative integers indexing 
the corresponding eigestate. Let us now discuss how to 
construct a VR starting from the considered linear limit. 
Intriguingly, utilizing the results of Ref. [3l| , this is pos¬ 
sible in the following way: in that work it was found 
that for the anisotropic trap with = 1 and ojz = 2/fc 
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the energy of a second radial excited state {m,n,k') with 
m + n = 2 and k' = 0 (i.e. E = 3 + 1/k) coincides with 
that of the fc-th excited state along the z-axis (0, 0, k) 
(also equal to 1 + (2/fc)(fc + 1/2) = 3 + 1/k). For such 
anisotropic traps one can construct stationary states with 
/c-parallel-VRs 


u{x,y,z) oc 


| 2 , 0 , 0 )+ 10 , 2 , 0 ) 

V2 


i\0,0,k), 


( 4 ) 


even at the linear limit. The RDS nature of the real part 
combined with the k oscillations of the imaginary part 
along the z-direction constitutes parallel VRs with alter¬ 
nating vorticity. Remarkably, this is a very natural 3D 
generalization of the 2D setting considered in Ref. [s^. 
In the above expression of Eq. (SD, we have used the no¬ 
tation 


|n,m, fe) = 

X eyzp[-{uJr[x^ + y'^) + ujIz^)/2], (5) 

corresponding to the eigenmode of energy En^m.k of the 
QHO where El represents the Hermite polynomials. 

Generalizing the approach of Ref. [341 ] by consider¬ 
ing anisotropic trap strengths and (including the 
isotropic one of Wi- = Wz), the energies of the RDS and 
the fc-th solitonic state will, respectively be, 

Erds = 3a;r H—and Egoi = + (fc + 

Hence, for kujz < 2ujr, the state of lower energy will be 
the fc-th soliton from which the fc-th VR will bifurcate, 
while for kujz > 2wr, the lower energy will belong to the 
RDS, with the bifurcation occurring from that state. 


1. Single vortex ring, ujz < 2uJr 


Following the bifurcation phenomenology of Ref. 

(see also Ref. (l^), we first consider the limit ujz < 2uJr 
which encompasses the isotropic trap, and proceed as 
follows. Take the fc-th soliton along the z-axis, focusing 
for now on fc = 1, and assign this as mode ui with energy 
El. Similarly, assign the RDS as mode U 2 with energy 
E 2 . If a novel state (in this case, the single VR) bifurcates 
from these two states when their phase difference is tt/2, 
then a general two-mode analysis [s^ predicts that the 
number of atoms N = J \u\'^dxdydz at the bifurcation 
critical point will be given by 




El — E 2 
hi — III ’ 


( 6 ) 


where for the above special case of the single VR we 
have El — E 2 = uJz — 2uJr, while In = f \ui\'^dxdydz and 
hi = J \ui\'^\u 2 \'^dxdydz are overlap integrals. Specif¬ 
ically, for the VR, the dark soliton state from which 
the bifurcation arises is Ui = |0,0,1), while the RDS 


IS U2 = 


| 2 , 0 , 0 )-|-| 0 , 2 , 0 ) 
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Furthermore, the general two-mode theory also pre¬ 
dicts the chemical potential value at which the bifurca¬ 
tion will occur, namely 


Mcr — El + IiiNcr, (7) 

where N^r is given by Eq. d®]). 

Interestingly, the above integrals can be computed an¬ 
alytically. In particular, for the case of the single VR 
with fc = 1, we find that 


hi 


3ujj.,yujz _ 
8y27r3/2 “ 12 


( 8 ) 


As a result, the explicit prediction for the bifurcation of 
the VR is that: 


Ar(v+i) 

'cr 

4-\/27r3/2(2a;7. — w^) 

(9) 

LVry/iOz 

.(VR,1) 

Mcr 

3 

II 

(10) 


This provides us with an explicit prediction for the bi¬ 
furcation of the first VR. 


2. Single vortex ring, uiz > 2uJr 


For completeness, we now consider the case of uiz > 
2ujr, where the VR now bifurcates from the RDS and the 
subscripts of Eq. ([51) must be exchanged 102. For this 
case, E 2 — E 1 = 2ujr—uJz and I 22 = / \u 2 \‘^dxdydz = 2 / 12 , 
which leads to the bifurcation point: 


Ar(VR.2) 

cr 




(VR,2) 

cr 


8-\/27r3/2(a;2 — 2u}r) 

UJr ^ UJ z 

5 

-UJz - LOr- 


( 11 ) 

( 12 ) 


The superscripts in Eqs. and CHD versus those of 
Eqs. (fTT|l and (fT^ are used to illustrate which state the 
VR bifurcates from, 1 is for the dark soliton and 2 rep¬ 
resents the RDS. 


3. Double vortex rings, tOz < o 


Similarly, for the double-vortex-ring (2VR) state Ei — 
E 2 = 2{uiz —ujr), while the overlap integrals are given by 


4l0Jr^/uh 

128y27r3/2 



(13) 


Hence, the corresponding prediction for the bifurcation 
point yields 


Ar(2VR.i) 

'^cr 

,,(2VR.1) 

A*cr 


256y27r3/2(w^ — OJz) 
290Jr\/‘jIz 
222^1^7. — 49 iOz 
58 ■ 


(14) 

(15) 


These predictions are valid for lOz < oJr- In this case, the 
2VR will bifurcate from the two-dark-soliton state of the 
form ui = jo,0,2), while the higher energy state in the 
two-mode analysis will again be the RDS. 
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4- Double vortex rings, lJz > uir 


(|2,0,0)-|0,2,0))/V2, one finds the bifurcation point: 


On the other hand, for tOz > oJr the prediction needs to 
be suitably modified. Since the 2VR now bifurcates from 
the RDS we again exchange the subscripts in Eq. (El- 
Using I 22 = (41w^^)/(128y27r3/2) (for U 2 = |0,0,2)) 
we finally find that: 


at{2VR,2) 

''cr 

(2VR.2) 

P-cr 


64-\/27r^/^ — LOr) 

BuJry'^ 

^Tujz 

lo ■ 


(16) 

(17) 


/V(2SV,1) 


p 


(2SV,1) 

cr 


4-\/27r^/^(2a;r — uiz) 
AuJr, 


( 20 ) 

( 21 ) 


for ujz < ‘i.ujr ■ The reverse trap anisotropies can similarly 
be explored. 


6. Comparison to other work 


An interesting feature is the following. For uJz < 2wr, 
the single VR bifurcates from the one-dark soliton (fc = 
1) and for uiz < uir the 2VR bifurcates from the two- 
dark-soliton state (fc = 2). In the “intermediate” case, 
ujr < oJz < 2uJr, the 2VR state bifurcates from the RDS, 
while the single VR bifurcates from the one-dark soli¬ 
ton. Finally, for 20;^ < both the single VR and the 
2VR states bifurcate from the RDS. Hence, as expected, 
^gVR. 2 ) ^ ^(yR. 2 )^ lower-energy single VR 

bifurcates, followed by the 2VR state. 

We remark in passing that, in line with the special-case 
results of Ref. , the single VR bifurcates from the lin¬ 
ear limit {N —>■ 0) in the anisotropic case of tUz = 2wr, 
as per Eq. (jO]), while the 2VR bifurcates from the lin¬ 
ear limit, precisely for the isotropic case of lUz = w,., as 
per Eq. m- One can similarly generalize these types 
of bifurcation considerations for all higher-order rings, 
providing an explicit set of predictions for their emer¬ 
gence in the vicinity of the linear limit. This bifurcation 
and stability analysis naturally explains why the decay 
of dark solitonic states yields the corresponding k-VR 
states, not only in numerical simulations [33|, but also in 
experiments [12| . 


5. Solitonic vortices 


Finally, and before proceeding with the stability anal¬ 
yses of the dark soliton and the resulting VR states that 
are the principal focus of the present work, let us compare 
our work to Ref. [s^. Our existence results bear signifi¬ 
cant resemblance to those of Ref. , although there are 
also nontrivial differences. For instance, in Ref. [s^ the 
focus was on bifurcations from the dark soliton (called the 
kink therein) state, whereas we especially focus on bifur¬ 
cations from various states that yield VR-like solutions 
(single and double VRs). In Ref. [s^, the z-direction 
is presumed to be homogeneous (i.e., untrapped), while 
here we consider the presence of a trap. For this reason, 
multi-VR states (such as those explored above) are not 
possible in the framework of Ref. [^. In fact, the 2VR 
state in that work is one in which both rings are in the 
same plane, while here we explore multiple non-coplanar 
rings in a stationary state, rendered possible by the addi¬ 
tional trapping along the z-direction. On a technical level 
too, the methodology of the reduction to a quasi-linear 
equation on the plane with an effective potential used in 
Ref. im, while ingenious, differs substantially from the 
two-mode approach utilized herein. Finally, and perhaps 
most importantly, a principal focus that stems from our 
existence findings is that, we also study the stability of 
the obtained solutions (see details in the following sec¬ 
tion). In Ref. [13, on the other hand, this is deferred to 
future studies. 


While the emphasis here is on VRs, partly to illus¬ 
trate that the different states obtained in Ref. [13 can 
be identified with the techniques presented herein, we 
also provide other special cases, namely the single- and 
double-solitonic vortex states. 

For the case uJz < uJr, following the same technique as 
above but with U 2 = |1,0, 0), we obtain: 


Ar(isv,i) 




(ISV.l) 

cr 


4-\/27r^/2((^^ _ 
OJr y/ ^ z 

5uJr 


(18) 

(19) 


B. Stability Analysis Near the Linear Limit 


We now move to the stability analysis in the vicinity of 
the linear limit. Using a Taylor expansion of the solution 
u = y/luQ + + ... and ^ = /ro + e/ri -I-... in Eq. ([T|) , 

where (p-o^uo) correspond, respectively, to the eigenvalue 
and eigenfunction of a state at the linear limit, we find 
at 0(e) the solvability condition: 


Ml 


\uo\'^dxdydz. 


( 22 ) 


for the single solitonic vortex (bifurcating from the dark 
soliton Ml = |0,0,1)). For the double solitonic vor¬ 
tex, again bifurcating from mi = |0 ,0,1) but with U 2 = 


This, in turn, allows us to specify e = (/r — ^o)/Mi (to 
leading order). Then, the spectral stability, as discussed 
in Ref. [13 , amounts to solving the BdG eigenvalue prob- 
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lem ("Ho + Mi)U = uU, with 


U 


'Ho 



where C = —(1/2)V^ + V{r) — ^o, while 


Hi 


2|moP-Mi “o 

-(wq)* /ri-2|uoP 


(23) 

(24) 


(25) 


where the star denotes complex conjugate. Here, we de¬ 
note oj as the eigenfrequency of a given eigenstate U with 
eigenvalue X = ito. The presence of a nonzero imaginary 
part of w or, equivalently, of a nonzero real part of A in 
our Hamiltonian system denotes the presence of a dy¬ 
namical instability. 


1. The dark soliton 

From the above formulation, it is straightforward to 
analyze the stability in the case of e —>■ 0. There, it is 
evident that the linearization spectrum consists of the 
diagonal contributions of Ho which consist of the spec¬ 
trum of L and of its opposite. Up to now, we have kept 
our exposition as general as possible, but from here on, 
we will focus on a specific state in order to showcase the 
relevant ideas more concretely. As our workhorse, we will 
use the dark soliton state uq = |0, 0,1), in a spirit similar 
to that of Ref. [s^ . but with a particular view towards 
the bifurcation of the VR state (as well as others such 
as the 2VR and the single and double solitonic vortex). 
Since C consists of the quantum harmonic oscillator, the 
spectrum for a dark soliton along the z-axis, when sub¬ 
tracting fiQ = + i^zl‘2., will be (see Eq. ([3])) 

u! = uJr {ri + m) + u}z{k — 1). (26) 

Recalling that n, m, and k are arbitrary non-negative 
integers, we now proceed to explore the relevant states 
and their multiplicities. The mode with n = m = k = 0 
corresponds to the negative energy im, or negative 
Krein signature [i^, mode. Such a mode, when reso¬ 
nant with another positive energy (or Krein signature) 
one, will give rise to complex eigenvalue quartets, while 
collisions of modes of the same energy will be inconse¬ 
quential towards changing the stability properties of the 
solution. The dark soliton state (with vanishing den¬ 
sity in the z = 0 plane) has a neutral mode, corre¬ 
sponding to {n,m,k) = (0,0,1), which has a wavefunc- 
tion identical to the solution itself. This mode corre¬ 
sponds to the phase or gauge (U(l)) invariance of Eq. ([T|). 
There are three additional Kohn modes corresponding 
to symmetry which are also left invariant. Specifically, 
(n, m, fc) = (0,0,2) corresponds to dipolar oscillations 
along the z-axis with frequency oj = while the modes 
{n,m,k) = (1,0,1) and {n,m,k) = (0,1,1) pertain to 


dipolar oscillations along the {x, y)-plane with frequency 
uj = UJr- To complete our discussion of modes with 
n -I- m -I- < 2, we need to account for five more modes. 

There are two degenerate modes (due to the radial invari¬ 
ance of the trap in the plane) of (n, m, k) = (1,0, 0) and 
(n, m, k) = (0,1, 0), with frequency a; = uor—uoz and three 
degenerate ones, (n, m, k) = (2,0, 0), (n, m, k) = (0, 2, 0) 
and {n,m,k) = (1,1,0), all of which have frequency 
Lj = 2ujr — uJz in the linear limit. 


2. The dark soliton, u)r = uoz 


The key question that subsequently arises is that of 
the fate of the eigenvalues described above, as e becomes 
finite, i.e., as we depart from the linear limit. To follow 
these eigenvalues, and given that the different degenera¬ 
cies also hinge on the specific values of uJr and we 
will use as our benchmark case the isotropic scenario of 
Wr = Wj = 1, where the choice of unity is made without 
loss of generality. In this case, the zero eigenvalue has a 
multiplicity of three. The gauge-invariance mode must 
remain at zero, and so too must the two modes (1,0,0) 
and (0,1,0) due to the spherical invariance of the dark 
soliton solution, but only when the trap is isotropic. 

We now move to the consideration of the principal 7 
modes (recall that this is really 7 pairs) at w = ±1. Out 
of these, the 3 dipolar modes will remain invariant. Then, 
however, 4 modes are subject to deviations, as soon as 
we depart from the linear limit. It turns out that the 
insightful work of Ref. [s^ has already computed one 
sub-manifold associated with such a bifurcation. In that 
work, the authors recognized that the eigenvector associ¬ 
ated with the RDS Ui = (0)'^ [see Eq. (1331) ] 

and the anomalous mode U 2 = (0, |0, 0,0))^ become res¬ 
onant, and identified the deviations by using degenerate 
perturbation theory, i.e., constructing the matrix M. with 


M,j = {U,\'Hi\Uj), (27) 


and identifying its eigenvalues for i,j = 1,2 for the above 
eigenvectors. That calculation, adapted to the present 
setting, can be rewritten as 


(3±^ 

12 



(28) 


(cf. Eq. (20) in Ref. [H). It is important to note that the 
near-linear prediction is that the dark soliton should be 
immediately unstable due to this resonant interaction, via 
an oscillatory instability and a quartet of corresponding 
eigenvalues. This is a general feature that dark soliton 
(and multi-soliton) states possess near the linear limit 
due to the degeneracy of their anomalous modes; cf. the 
work of Ref. [^. However, in that work these instability 
“bubbles” were shown to terminate at some finite value 
of /r, hence it is of interest to explore whether a similar 
feature arises here, a question that we will address below 
numerically. 
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On the other hand, in the work of Ref. [s^, an addi¬ 
tional (and, as we will see, important) sub-manifold of 
eigenvectors was not considered, namely that of C/i = 
( |2.0.0M0.2,0) ^q^t Qf ^ (|1,1,0),0)^. For this 

subspace, we find that the corresponding deviation of the 
eigenfrequency from the linear limit —again, obtained 
via the degenerate perturbation theory of Eq. dill)— 
yields: 

Given the decreasing trend of both of these eigenvalues, 
it is natural to expect that at some finite value of /r, they 
will hit the origin of the spectral plane. As a preamble 
to our numerical computations of the next section, it is 
then relevant to consider what the outcome of such a col¬ 
lision will be. Connecting these findings with our bifurca¬ 
tion theory results of the previous subsection, we appre¬ 
ciate that these collisions should be, in fact, what leads 
to the corresponding bifurcations and destabilizations of 
the dark soliton along the z = 0 plane. In particular, the 
mode Ui = (0)'^ effectively corresponds to 
the RDS. Its “collision” with the origin and subsequent 
destabilization of the dark soliton along the z = 0 plane 
suggests that beyond this threshold the planar dark soli¬ 
ton and RDS mesh, which, as we have discussed before, 
produces the single VR state. In the same spirit, we can 
see that the collision of Gi = ^ q)^ with the 

origin will produce a meshing with the planar dark soli¬ 
ton and lead to the bifurcation of the two solitonic vortex 
state, discussed in the previous subsection. In fact, as¬ 
suming the prediction of Eq. (1291) to be useful beyond 
its realm of validity (and until w = 0) yields a critical 
point for the relevant bifurcation at ^ = 4, which coin¬ 
cides with the prediction of p,GSV,i) ([2T]) . We will 

numerically examine the validity of these predictions in 
the next section. 


3. The dark soliton, ujr = 2lOz 

In principle, the approach adopted above (and the as¬ 
sociated eigenvalue count and deviations from the linear 
limit) can be used for any state bifurcating from the lin¬ 
ear limit. However, obviously, the more complicated the 
original state, the more difficult it becomes to account 
for all the relevant eigenvalues. Since our focus here is 
on VRs and their emergence, we will also give an ad¬ 
ditional example of the case of utr = 2 (jJz = 1 (more 
generally bearing in mind the case oi LOr ^ toz), which 
will be considered in our numerical eigenvalue analysis 
below. Once again, in this case, 4 modes remain invari¬ 
ant, namely one at w = 0 (phase invariance), and at 
u) = LOz and (double) to = uor due to the dipolar oscil¬ 
lations. The manifold of the modes (1,0,0) and (0,1,0) 
with frequency lo = uOr—tOz sA, the linear limit leads to the 


prediction (using again degenerate perturbation theory) 


UJ = (lOr — LOz) — 


1 


fi- {LOr + -LOz) 


(30) 


On the other hand, the anomalous mode of indices 
(0,0,0) is theoretically predicted to move according to 


LO = iOz 


1 

6 


^1- {LJr + -LOz) 


( 31 ) 


Then, the sub-manifold of eigenvectors with Ui = 
( | 2 .o,oMo, 2 ,o) ^ 2 ^^ 0), 0)^ produces two 

coincident eigenfrequencies (actually two pairs) with 


LO = {2lOr — LOz) — 


^Ji- {LOr + -LOz) 


(32) 


However, we can see that there are additional eigen¬ 
values, in this case, that acquire comparable values and 
therefore manifolds of larger sum n -|- m -I- fc need to be 
considered (up to now we had restricted considerations to 
n+m+k = 2). Among the eigenvalues with n-|-m-|-fe = 3, 
the case of (0, 0, 3) will have an eigenfrequency of w = 1 
in the present setting. The perturbative calculation in 
this case yields a prediction of 


LO = 2lOz 


1 

12 


fi- {LJr + -LJz) 


(33) 


Even more complicated are the degeneracies occurring 
at a; = 1.5. In addition to Ui = g)"^ 

and of U 2 = (|1,1,0),0)^ considered above, there is 
the double degeneracy of C/i = g)"^ and 

U 2 = (|0,0,4), 0)"^ and finally that of Ui = (|1,0, 2), 0)^ 
and U 2 = (|0,1, 2), 0)"^. The latter, leads to the degener¬ 
ate eigenfrequencies (again, two pairs) of the form: 


LO = {Wr + LOz) 


5 

12 


(LJr+ -LOz) 


(34) 


The former eigenvalues (i.e., the RDS and the one with 
fc = 4) are only resonant when = 2lOz. In this partic¬ 
ular case, we can obtain the corresponding eigenfrequen¬ 
cies 


LO 


1.5- 


53±7v^ 

192 



(35) 


We believe it is clear that while the relevant methodol¬ 
ogy is entirely general, the logistics of its application vary 
from case to case and can be fairly complex. Eor this 
reason, we now corroborate and complement our analy¬ 
sis with detailed numerical computations in the follow¬ 
ing section. The above two case studies, of the isotropic 
regime and of lj^ = 2lOz, will operate as our benchmarks. 
We will also numerically explore other settings such as 
LOz = ‘^LOr as well as, importantly, the spectra of our 
states of particular focus, namely the bifurcating (single 
and double) VRs. 
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FIG. 1: (color online) Spectrum for the dark soliton station¬ 
ary state for the isotropic case, ujr = uJz. Depicted are the 
eigenfrequencies stemming from the BdG analysis, as a func¬ 
tion of the chemical potential. Real parts of the eigenfrequen¬ 
cies are denoted by black cross points, while the imaginary 
parts of the eigenfrequencies are denoted by red plus points. 
The blue solid lines correspond to theoretical predictions in 
the small-interaction, linear, limit (cf. text). 


III. NUMERICAL RESULTS 

We identify the stationary solutions by solving 
the Gross-Pitaevskii equation using a Newton-Krylov 
scheme [1^. The Bogoliubov linearization spectrum is 
then obtained by utilizing a spectral basis of noninteract¬ 
ing modes; at least 800 are used for each result herein. 
From a numerical standpoint, the problem is reduced to 
two dimensions by utilizing a Fourier-Hankel method 
(see also Ref. [4^), which is made possible by the az¬ 
imuthal symmetry of the trap in the {x, y) plane. 

Recall that in Bogoliubov theory every eigenfrequency 
belongs to a pair, having solutions of opposite sign. From 
here on, we simplify our discussion by only plotting half 
of the spectrum, one eigenfrequency from each of these 
pairs. 


A. The dark soliton, Ur = 

We begin our discussion by considering dark soliton 
stationary state for the isotropic case, uJr = oJz = 1. 
The corresponding results are depicted in Fig. [T] a typ¬ 
ical example of the solution itself (both planar cuts, as 
well as a full 3D density isocontour plot) is shown in 
Fig.H A first observation is that our analytical predic¬ 
tions seem to agree well in the linear limit with the corre¬ 
sponding numerical results. More specifically, we observe 
that the complex quartet of Eq. (1^ indeed destabilizes 



FIG. 2: (color online) Dark soliton state for the isotropic case 
of ojr = Wz = 1 and y = 12. The top four panels show the 
modulus of the solution (left subpanels) and their respective 
argument 6 (right subpanels). The top row illustrates the 2 = 
0 plane, where the density vanishes, and the middle row shows 
the X = 0 plane. The bottom panel shows two isocontour 
density plots for this state (the contours correspond to iso¬ 
density surfaces at the maximum density divided by 2.5 and 
1.5). 


the planar dark soliton in the linear limit, as was origi¬ 
nally observed in Ref. [3^ . However, as we depart from 
this limit, the phenomenology reported in Ref. ap¬ 
pears to arise, namely the resonant interaction of the 
RDS mode and of the anomalous mode ceases. There¬ 
after, the anomalous mode maintains a roughly constant 
frequency, while the RDS rapidly decreases in frequency, 
eventually colliding with the origin (i.e., with a; = 0) 
for a value oi y = 4.05 (very close to the theoretical 
predictions of ^ = 4, see Eq. m)- Indeed, beyond 
this critical point, this eigenfrequency becomes imagi¬ 
nary (equivalently, the eigenvalue becomes real), giving 
rise to a symmetry-breaking pitchfork bifurcation. The 
daughter branch emerging from this bifurcation (per our 
analysis), is the single VR (of charge either -|-1 or — 1). 
This will be further corroborated by the computation of 
this state below. 

Eurthermore, we can observe a good agreement also 
for the two degenerate eigenfrequencies, stemming from 

























FIG. 3: (color online) Spectrum for the dark soliton station¬ 
ary state for = Suj = 1. Points have the same meaning as 
in Fig. [U 


Ui = and U 2 = (|1,1, 0), 0)"^ (where U 

is defined in Eq. (l23l) . in accordance with Eq. (|2^ . In 
this case too, the decrease in frequency eventually leads 
to a zero-crossing and a destabilization of the dark soliton 
state in favor, in this case, of a double solitonic vortex 
(2SV) state. It is important to note that, contrary to the 
homogeneous (along z) case reported in Ref. [32|, here the 
bifurcation of the 2SV and that of the single VR do not 
occur at the same value of the chemical potential. A con¬ 
sequence of this is that the dark soliton has already been 
destabilized by the bifurcation of the 2SV state when the 
VR bifurcates, and hence the VR in this isotropic case 
is expected to inherit this weak instability close to that 
limit. Note that there are three modes that are invariant 
at w = 0, one of which is due to the gauge invariance, and 
two have frequency Wr ~ Wz = 0, reflecting the rotational 
invariance of the planar dark soliton. 


B. The dark soliton, ujr = 2a;z 

We now turn to the investigation of our second analyt¬ 
ically examined case, namely of uJr = ‘2-uJz = 1. The 
comparison between the analytical prediction and the 
numerical results for the stability spectrum is depicted 
in Fig. [3] Remarkably, we can see in this case too — 
despite the considerable complexity of the bifurcation 
diagram— the high accuracy of our eigenvalue predic¬ 
tions near the linear limit. While in this case there is 
no quartet from the linear limit, the instability emerges 
due to the rapid decrease of the modes (accurately pre¬ 
dicted by Eq. dsnu, near the linear limit) pertaining to 
the manifold of (1, 0,0) and (0,1, 0). We remark that, in 
the isotropic case, this pair of modes was equi-energetic 



FIG. 4: (color online) Spectrum for the dark soliton station¬ 
ary state for ujr = 0 ^ 2/2 = 1. Otherwise, similar to Fig. [T] 
Here the ring dark soliton and the planar dark soliton have 
the same energy, and hence the VR can bifurcate immediately 
from the linear limit. The principal instabilities of the dark 
soliton in this case arise due to higher modes (with n-l-m = 3), 
although lower-order modes contribute to oscillatory instabil¬ 
ities associated with complex eigenfrequency quartets. 


with the dark soliton, allowing the construction from the 
linear limit of the single solitonic vortex. Here, however, 
it is at ^ = 2.58 that the relevant instability emerges, 
giving rise to the solitonic vortex state. Notably, this 
also agrees with the analytical prediction of /r = 2.5, 
see Eq. (HU). The only other instability that can be ob¬ 
served is given by the collision of the RDS, predicted by 
Eq. (l35l) . with the anomalous mode of Eq. (1311) . This, 
once again, leads to a quartet of eigenfrequencies and 
only for considerably larger values of ^ (not included in 
the hgure) will the VR bifurcate. Additional predictions 
such as the manifold of Ui = Q)^ and of 

U 2 = (|1,1, 0), 0)^ through Eq. (I33|) . as well as the higher 
index (n+m+k) of Eqs. ( 1331 ) and ( 1331 ) are also reasonably 
accurately captured. 


C. The dark soliton, ujz = ‘2ujr 

Finally, we also touch upon the case of uJz = = 2 

—cf. Fig. 01 In this case, the situation is considerably 
more complex and numerous instabilities appear to arise. 
Nevertheless, in this case too, our understanding of the 
linear limit may provide a reasonable set of guidelines 
for understanding the relevant phenomenology. In this 
case, we observe four instabilities associated with imagi¬ 
nary eigenfrequencies and three associated with eigenfre¬ 
quency quartets. Among the latter, the first to arise is 
the degenerate pair (the cross states) of eigenfrequencies, 
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M 

FIG. 5: (color online) Single VR for uir = oJz- Similar to 
that of Fig. [T] the present figure shows the BdG spectrum 
(this spectrum coincides with that of the dark soliton when 
the VR bifurcates from the it at /r = 4.05; cf. Fig. [!}. Notice 
that despite oscillatory and imaginary eigenfrequency-related 
instabilities (discussed in the text) for small values of fi, the 
VR is stabilized in the large chemical potential limit. 


emerging from the origin with negative energy and collid¬ 
ing with a higher-order (in n + m + k) mode at ^ = 5.59. 
The second quartet begins at /r = 6.21 and involves the 
RDS (arising from the origin with negative energy) col¬ 
liding with a higher energy mode. It should be remarked 
that the presence of the RDS at the linear limit enables 
(as discussed earlier) the bifurcation of the VR already 
at the linear limit. This, in turn, implies that the VR 
may be expected to be robust in this case (however, see 
the detailed results below). Finally, the third collision in¬ 
volves the degenerate modes (1,0,0) and (0,1,0) which 
are also now anomalous due to their negative energy and 
are growing from the linear limit. This pair collides with 
higher order modes and a quartet develops for fj. > 6.69. 
The first two imaginary instabilities arising are, surpris¬ 
ingly, not caused by lower-order modes, but rather from 
the higher-order ones (3, 0,0) and (0,3, 0), as well as from 
(1, 2,0) and (2,1,0). Both of these sets of eigenfrequency 
pairs decrease, as we depart from the linear limit, leading 
to an imaginary eigenfrequency instabilities for ^ > 5.24 
and fj, > 5.44, respectively. Hence, we see that in this 
case too, we can form a qualitative picture of the stability 
landscape merely by knowing and suitably appreciating 
the linear eigenfrequency/eigenmode picture. 


D. The single vortex ring, cur = oJz 

We now turn our attention to the case of the VR, first 
considering the isotropic trap; see Fig. [5] for the BdG 
spectrum and Fig. [5] for an illustration of the relevant 



FIG. 6: (color online) Single VR state for ujr = uJz = 1- Sim¬ 
ilar layout as in Fig. [2l The cross sections of the top two 
rows and the density isocontour of the middle panel show a 
prototypical example of the bifurcating branch (for /r > 4.05) 
of the single VR in the large chemical potential limit (here at 
^ = 12), where the relevant state has been stabilized against 
small-amplitude perturbations. The darker (green) isocon¬ 
tour corresponds to an iso-density at the maximum density 
divided by 8 plotted inside the bulk of the cloud. The bottom 
two rows correspond to a VR solution right after its bifur¬ 
cation at /r = 4.05, clearly illustrating its ring dark soliton 
character in the {x, y) plane and its “emerging vorticity” in 
the phase which is reminiscent of a vortex dipole (despite the 
resemblance with a dark soliton in the density). 
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state far from (top two rows) and near to (bottom two 
rows) its bifurcation point. As expected, the eigenfre- 
quency pertaining to the RDS Ui = ^ q)^, 

grows along the real axis. Recall that, for the dark-soliton 
stationary state, the RDS was the mode that became un¬ 
stable and gave rise to the VR. However, for larger values 
of the chemical potential (i.e., for fi G [4.35, 5.14]), and in 
a way reminiscent to what happens for the dark soliton 
in the vicinity of the linear limit, the RDS collides with 
the anomalous, negative energy mode U 2 = (0, |0,0,0))"^, 
giving rise to an oscillatory instability (and an associated 
eigenfrequency quartet). This instability disappears for 
larger values of ^ > 5.14. 

Importantly, since the VR bifurcated from the dark 
soliton state (both with ujr = ojz) then at the point of 
bifurcation their spectra should be identical, cf. Figs. [1] 
and E] at /r = 4.05. Recall, though, that at the bifur¬ 
cation point another instability was already present in 
the spectrum of the dark soliton due to the prior bi¬ 
furcation of the 2SV. Consequently, the VR is endowed 
with this instability from birth. However, upon increase 
of the chemical potential, we observe that the eigenfre¬ 
quency associated with this inherited instability (cor¬ 
responding to a doubly degenerate subspace associated 
with Ui = ( |2^FOMo.2,o) ^Q^r ^ (|1,1,0),0)^, 

i.e., the “cross states”) decreases and eventually leads to 
a complete stabilization of the VR for /i > 5.54. It is due 
to this stabilization that such VR patterns can be ob¬ 
served robustly in the large chemical potential (Thomas- 
Fermi) limit. In the latter limit, the VR acquires particle¬ 
like characteristics; see, e.g.. Ref. fo r the dynamics 
of a single VR inside a trap and Ref. for the interac¬ 
tions between multiple VRs. The combination of these 
features and the connection of the spectral features of 
the VR with the above particle-like characteristics will 
be explored separately in a future publication. 


E. The single vortex ring, ljz = 2(jJr 

We now briefly examine the VR for the special case 
where it bifurcates from the linear limit, namely when 
LOz = ‘2-uJr (and hence the energies of the dark soliton 
and the RDS become degenerate, enabling the construc¬ 
tion of a VR already at the linear limit, in accordance 
with Eq. @ ). The corresponding BdG spectrum is de¬ 
picted in Fig. [T] Interestingly, here, we observe that 
the VR is, in fact, unstable immediately (i.e., as soon 
as we depart from the linear limit). This instability 
is manifested through a degenerate imaginary pair of 
eigenfrequencies, again associated with the cross states 
Ui = of ^ (|1,I,0),0)^. In 

this limit, the latter states also bifurcate from the lin¬ 
ear limit and apparently have lower energy than the 
VR. We note, however, that as we proceed to larger val¬ 
ues of the chemical potential, the instability appears to 
(weakly) decrease and, hence, VRs may again be long- 



FIG. 7: (color online) Spectrum for the single VR when = 
LtJzl2 = 1, similar to Fig. [5] Here, the VR emerges from the 
linear limit. 





FIG. 8: (color online) Spectrum for the double vortex ring 
(2VR) for the isotropic case Ur = oJz = 1. Points have the 
same meaning as previous figures. 


lived in the Thomas-Fermi limit of large /r. Addition¬ 
ally, it is worthwhile to note that once again the colli¬ 
sion of Ui = ( l^’0’°)Jdo,2,o) ^o)r anomalous mode 

U 2 = (0, jo, 0,0))^ leads to a resonance and an oscillatory 
instability for the interval 5.78 < /r < 8.00. 

We also remark that, for this stationary state as well, 
the relevant eigenvalue count can be performed. Consider 
the eigenvalues of the operator £ — /tq (with /tq = 40;^.) 
to he uJr{n + m + 2k) — 2ujr- We find that there are (a) 
4 modes at the spectral plane origin (3 climb from the 
origin with increasing as discussed above, and one re- 
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FIG. 9: (color online) Double VR (2VR) state for the isotropic 
case (jJr = (jJz = 1 and fj. — 12. Similar layout as in Figs. [2] and 
[6l This state bifurcates from the linear limit at /r = 7a;r/2. 


mains at the origin due to gauge invariance), there are 
(b) 8 modes at w = uJr, two of which are anomalous, 
while there are (c) 10 modes at uj = 2ujr, one of which is 
anomalous and gives rise to the oscillatory instability dis¬ 
cussed above. Hence, given the technical complications of 
considering the relevant perturbative analysis (and also 
the qualitative understanding afforded to us by means of 
the above analysis), we will not pursue this further here. 


F. The double vortex ring, = oJr 

Finally, we also examine the double-vortex ring (2VR) 
state, bifurcating from the linear limit of fiQ = luJrl‘2- in 
the isotropic case of see Figs. [5] and 111 Here 

the situation is even more complicated, with 6 degener¬ 
ate eigenfrequencies at w = 0, 13 modes (3 anomalous 
ones) at O’ = uir, and so on. Hence, we will again re¬ 
strict our considerations to some qualitative remarks. We 
note that the state will be immediately unstable in this 
isotropic limit, due to an imaginary eigenfrequency pair 
(the cross states), bifurcating from 0 as soon as the 2VR 
state emerges. Moreover, the resonance of the degener¬ 
ate anomalous pair at uj = uir with the corresponding 
positive (same) energy modes leads to a scenario again 


reminiscent of Ref. in that an oscillatory instabil¬ 
ity emerges from the linear limit, although it is termi¬ 
nated at fi = 4.47. An additional collision of a mode 
bifurcating from the origin with the anomalous mode 
U 2 = (0, |0,0, 0))^ leads to an additional eigenfrequency 
quartet for jj, > 4.96. A key observation here is that these 
instabilities (both the imaginary one and the oscillatory 
one for /i > 4.96) were found to persist throughout the 
interval of parameters considered herein. However, for 
large values of yc, we observe a weak tendency for both of 
these instabilities towards decreasing growth rates. The 
latter feature appears to be qualitatively consonant with 
the observations of Ref. 0, suggesting that, for inter¬ 
mediate parameters, 2VR states were less robust than for 
large values of the corresponding nonlinearity-controlling 
parameter. Once again, the theme of the Thomas-Fermi 
limit will be the subject of separate work, focusing on the 
dynamics and interactions of VRs as particle entities. 


IV. CONCLUSIONS AND FUTURE 
CHALLENGES 


In conclusion, we have systematically examined the full 
eigenvalue spectrum by linearization around a dark soli- 
ton in a 3D setting, building on considerable earlier work. 


most notably that of Ref. [33| . Extending that work, we 


have shown how to predict not only the formation and 
bifurcation of the (single) vortex ring state, but also of 
other states, such as the double-vortex-ring, the single- 
and the double-solitonic vortices, among others. We have 
explained how the emergence of these states can be pre¬ 
dicted on the basis of a Galerkin-type, two-mode, ap¬ 
proach. This is analogous to how a vortex pair, and 
more complex states involving multiple vortices, were 
predicted to arise in 2D (see, e.g.. Refs, [sl, [s^), a 
step that subsequently led to their experimental real¬ 
ization [ 4 ^. We have explored the conditions necessary 
for the vortex ring to emerge through a bifurcation from 
the planar dark soliton state, and explained when it can 
arise from the linear limit (when the energy of the ring 
dark soliton and that of the dark soliton coincide). Con¬ 
versely, the vortex ring may arise through a bifurcation 
from the ring dark soliton, when the latter possesses lower 
energy than the planar dark soliton state. We have also 
highlighted similarities and differences from important 
related work such as the contribution of Ref. . 

In addition to giving a qualitative characterization of 
the relevant spectra and bifurcations, in the vicinity of 
the linear limit, we also employed perturbation theory 
(often in its more complex degenerate form) in order 
to quantitatively characterize the eigenvalues responsi¬ 
ble for the relevant instabilities and the emergence of 
new branches. 

Upon obtaining a systematic understanding of when 
the vortex rings arise, we turned our attention to their 
spectral stability properties, by numerically solving the 
Bogoliubov-de Gennes equations. We thus found that at 
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the isotropic limit, the single vortex ring may be unsta¬ 
ble “at birth”, but can be stabilized at higher chemical 
potentials. In the anisotropic case, where the vortex ring 
state emerges from the linear limit, it is also immediately 
unstable but this instability can be relatively weak in dif¬ 
ferent parametric regimes, such as that of large chemical 
potential. Similar features were found for the double vor¬ 
tex ring, in qualitative agreement with earlier dynamical 
observations [3l|. 

We believe that this study paves the way for a deeper 
understanding of such vortex ring states, offering an un¬ 
precedented view of their spectral features. In this light, 
there are numerous aspects worthy of further exploration. 
A more technical example involves the attempt to sys¬ 
tematically characterize the spectrum of multiple vortex 
rings in the special cases where they emerge in the linear 
limit, i.e., at = 2ujr/k for the fc-th vortex ring state. 
A more intriguing aspect for our considerations is to ex¬ 
plore the opposite limit more systematically, namely the 
Thomas-Fermi realm, where the single and multiple vor¬ 
tex rings possess particle characteristics. While in this 
area there is some work already done, it is rather incom¬ 
plete. For instance, in Ref. [^, a theoretical approxi¬ 
mation of the vortex ring internal mode spectrum was 
obtained at this “particle limit”, but it was never tested 
against the full Bogoliubov-de Gennes spectrum or the 
three-dimensional GPE. On the other hand, in Ref. (d^ . 
a particle picture is derived and favorably compared to 
the GPE results, but this is only done in the absence of 
a trap. Obtaining a conclusive spectral picture for single 
and multi-vortex rings in the presence of a trap would cer¬ 
tainly be of paramount importance for our understanding 


of the particle-like character of these complex states and 
of their interactions within trapped BECs. Finally, once 
these aspects are addressed, one can think of not only 
examining the role of thermal and/or quantum fluctua¬ 
tions on these rings, but also importantly of generaliz¬ 
ing them in multi-component systems, where elaborate 
spinorial states exist in the form of skyrmions (50j | and 
monopoles m, among many others. Such studies are 
currently in progress and will be reported in future pub¬ 
lications. 
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